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CONFERENCE  PAPER 

This  paper  studies  the  inference  of  space  object  mass,  which  is  made  possible  due  to  the  coupled  influence  of  solar 
radiation  pressure  (SRP)  acceleration  on  the  orbit  of  satellites  and  their  observed  brightness.  This  effect  takes  time 
to  be  observed  in  optical  angle  measurements  given  the  combination  of  a  priori  kinematic  state  uncertainties  and  the 
magnitude  of  this  effect  relative  to  them  and  the  sensor  data  noise.  Therefore,  multiple  nights  of  observations  are 
typically  required  to  extract  this  “weak”  signal  from  collected  measurements.  From  angles  data  alone,  only  effective 
albedo-area-to-mass  can  be  estimated  since  this  term  appears  in  the  SRP  acceleration  equation,  but  when 
photometric  data  is  fused  with  the  astrometric  angle  measurements,  it  provides  observability  of,  and  thus  constrains, 
the  albedo-area  estimates.  This  inferred  constraint  makes  mass  the  most  open  degree  of  freedom  and  thus  the  fused 
data  eventually  informs  the  filter  of  the  mass.  The  observability  of  albedo-area  products  is  provided  by  the 
photometric  brightness  measurements,  since  the  brightness  of  the  space  object  is  a  strong  function  of  the  albedo- 
areas.  However,  the  relationship  between  the  albedo-areas  and  both  the  photometric  return  and  SRP  involves 
knowledge  of  the  Bidirectional  Reflectance  Distribution  Function  (BRDF)  for  the  surface  of  the  space  object.  If  the 
BRDF  in  the  photometric  measurement  model  and  the  BRDF  in  the  SRP  model  are  not  consistent  with  each  other, 
then  the  resulting  estimated  albedo-areas  and  mass  are  inaccurate  and  biased.  This  work  studies  the  use  of 
physically  consistent  BRDF-SRP  models  for  mass  estimation.  Simulation  studies  are  used  to  provide  an  indication 
of  the  benefits  of  using  these  new  models.  An  unscented  Kalman  filter  approach  that  includes  BRDF  and  mass 
parameters  in  the  state  vector  is  used.  The  full  set  of  estimated  parameters  includes  position,  velocity,  attitude, 
angular  rates,  mass,  exponential  factor  (parameter  in  Ashikhmin-Shirley  BRDF  related  to  sharpness  of  specular 
reflection),  specular  coefficient  and  diffuse  coefficient.  The  challenge  of  adding  these  additional  parameters  is  the 
fact  that  they  are  constrained,  where  all  must  be  positive  and  the  specular  and  diffuse  coefficients  must  be  less  than 
unity.  Two  approaches  are  adopted  to  account  for  the  constraints;  the  first  approach  projects  the  sigma  points  onto 
the  constraint  boundary  if  any  violate  the  constraints  on  the  parameters,  and  the  second  approach  defines  proxy 
parameters  that  are  unconstrained  version  of  the  original  parameters.  The  results  for  estimating  mass  are  promising 
and  show  that  the  addition  of  consistent  BRDF-SRP  models  benefits  the  estimation  process. 

1.  INTRODUCTION 

Deep  space  optical  surveys  of  near  geosynchronous  (GEO)  objects  have  identified  a  class  of  high  area-to-mass  ratio 
(HAMR)  objects  [1].  The  exact  characteristics  of  these  objects  are  not  well  known  and  their  motion  pose  a  collision 
hazard  with  GEO  objects  due  to  the  solar  radiation  pressure  (SRP)  induced,  large  variations  of  inclination  and 
eccentricity.  These  objects  are  typically  non-resolved  and  difficult  to  track  due  to  dim  magnitude  and  dynamic 
mismodeling.  Therefore,  characterizing  the  large  population  of  HAMR  objects  in  GEO  orbits  is  required  to  allow 
for  a  better  understanding  of  their  origins,  and  the  current  and  future  threats  they  pose  to  the  active  space  object  (SO) 
population. 

For  establishing  unique  identification  and  discrimination  of  SOs,  past  researchers  have  searched  for  defining 
attributes  that  are  unique  to  a  given  SO.  These  attributes  can  then  be  used  to  study  classes  of  objects  or  at  the  very 
least  correlate  observations.  An  obvious  attribute  that  can  be  assigned  to  an  SO  is  its  orbital  state.  Reference  [2] 
studies  the  process  of  correlating  a  sequence  of  optical  measurements  to  a  given  SO  using  constraints  on  the  orbital 
states  derived  from  the  measurements.  Although  orbital  states  do  provide  an  attribute  to  identify  an  SO,  for  the  case 
of  HAMR  objects  and  objects  with  sparse  observations  it  may  be  difficult  to  identify  an  SO  with  just  this  attribute. 
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In  order  to  understand  the  nature  and  eventually  the  origin  of  these  objects,  physical  characteristics  such  size,  shape 
and  material  are  required.  These  characteristics  can  serve  as  information  on  the  shape  and  the  attitude  motion  of  a 
debris  piece.  They  can  be  obtained  from  existing  data  sources  such  as  photometric  and  astrometric  data  [3],  To 
investigate  the  surface  material  properties  reflectance  spectroscopy  data  can  be  gathered  [4].  Reflectance 
spectroscopy  or  multi-band  photometric  data  can  be  used  to  study  non-resolved  debris  objects  and  determine  their 
material  composition  [5],  The  measurements  made  from  debris  objects  can  be  compared  to  laboratory  derived 
observations  to  determine  material  composition  of  the  debris  [6-8],  The  photometric  return  depends  on  a  number  of 
parameters  such  as  shape,  size,  material  and  observation  geometry.  Therefore,  extracting  material  composition  is 
difficult.  Reference  [9]  studies  the  light  curves  of  different  uncontrolled  SOs.  From  this  study  it  is  clear  that  there 
are  unique  features  in  these  light  curves.  But  reflectance  parameters  and  albedo  estimates  may  also  provide  useful 
attributes.  Size  estimates  have  been  derived  in  the  past  using  radar  measurements  [10],  but  for  objects  in  higher 
orbits  these  types  of  measurements  may  not  be  available.  Reference  [11]  develops  a  technique  for  estimating  the 
intrinsic  but  rough  size  and  albedo  estimates  of  orbital  fragmentation  debris  from  optical  measurements.  These 
albedo  values  provide  a  very  useful  attribute,  but  the  estimation  technique  assumes  that  the  objects  are  spherical  in 
shape.  This  assumption  limits  the  accuracy  of  this  method. 

Features  may  also  exist  in  the  light  curve  measurements.  For  example  some  objects  have  been  identified  by  the 
shape  and  features  in  their  light  curves.  Reference  [12]  shows  how  glint  data  can  be  used  to  determine  the  solar 
panel  off-set  and  improve  the  detectability  of  SOs.  Glint  data  however  is  very  much  a  function  of  the  geometry  and 
orientation  of  the  SO,  since  glint  events  only  occur  at  given  reflection  angles.  Therefore,  to  process  this  data  attitude 
and  shape  must  be  accounted  for  in  the  estimation  process.  Reference  [12]  does  not  use  estimation  techniques  to 
account  for  the  attitude  dependency,  but  rather  makes  qualitative  conclusions  from  glint  data.  Also  Ref.  [13]  uses 
optimal  control  methods  to  determine  reachability  and  correlate  observations;  knowledge  of  mass  can  greatly  assist 
in  these  types  of  calculations. 

Estimating  the  dynamic  characteristics  of  a  HAMR  object  using  light  curve  and  astrometric  data  can  allow  for  mass 
parameters  to  be  observable.  Estimating  mass  for  HAMR  objects  can  help  in  the  development  of  a  detailed 
understanding  of  the  origin  and  dynamics  of  these  objects.  It  has  been  shown  that  the  SRP  albedo  area-to-mass 
ratio,  CrA/m,  is  observable  from  angles  data  [14]  through  the  dynamic  mismodeling  of  SRP  forces.  Reference  [14] 
conducts  a  study  with  simulated  and  actual  data  to  quantify  the  error  in  the  estimates  of  CrA/m  and  good 
performance  is  found  using  data  spanning  over  a  number  of  months.  Also  Ref.  [15]  shows  that  orbital,  attitude  and 
shape  parameters  can  be  recovered  with  sufficient  accuracy  using  a  multiple -model  adaptive  estimation  approach 
coupled  with  an  Unscented  Kalman  Filter  (UKF).  This  approach  works  reasonably  well  but  requires  that  the  area- 
to-mass  ratio  be  known  a  priori.  The  purpose  of  this  work  is  to  show  that  since  CrA/m  is  observable  from  angles 
data  and  shape/albedo  properties  are  observable  from  photometric  data,  then  by  fusing  these  data  types  mass  can  be 
extracted  with  reasonable  accuracy. 

The  area-to-mass  ratio  of  a  space  object  can  also  provide  an  attribute.  This  ratio  is  very  important  for  predicting  the 
orbital  motion  of  the  SO  under  the  influence  of  SRP.  Since  SRP  is  the  primary  nonconservative  disturbance  force  for 
GEO  objects  it  is  a  very  important  dynamical  parameter  for  this  class  of  objects.  If  the  albedo  of  the  object  is  known 
it  may  be  possible  to  determine  the  mass  from  this  parameter.  A  similar  approach  has  been  applied  to  low-Earth 
orbiting  objects  using  the  drag  ballistic  coefficient  [16],  Knowing  the  mass  of  a  space  object  can  be  very  valuable. 
For  example  one  can  examine  breakup  fragments  by  determining  the  masses  of  each  fragment  [17],  This  work  will 
look  at  methods  for  determining  the  mass  of  an  SO  from  photometric  and  astrometric  data  sources. 

The  organization  of  this  paper  is  as  follows.  First,  the  methods  used  to  recover  mass  are  discussed.  Then  the 
models  used  for  SO  shape,  orbital  dynamics  and  attitude  dynamics  are  discussed.  Following  this  a  description  of  the 
measurement  models  used  in  this  paper  is  given.  Next,  a  review  of  the  constrained  UKF  approach  is  provided. 
Finally,  simulation  results  of  the  mass  and  albedo-area  estimation  approach  are  provided. 

2.  DETAILS  OF  APPROACH 

The  mass  estimation  approach  used  in  this  paper  depends  on  the  fusion  of  the  astrometric  and  photometric  data. 
Since  both  of  these  data  sources  have  an  albedo-area  dependency  then  mass  can  be  separated  from  albedo-area 
through  data  fusion.  For  example  consider  the  following  simplified  relationships  for  the  acceleration  and  observed 
flux  of  a  diffuse  sphere: 
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where  C,  is  the  reflectance  coefficient  of  the  sphere,  A  is  the  sphere’s  cross  sectional  area,  m  is  the  mass,  Fsun  is  the 
total  solar  flux  over  all  wavelengths,  c  is  the  speed  of  light,  Cdijf  is  the  diffuse  coefficient  of  the  surface,  Fsu no/jt 
represents  the  diffuse  BRDF  with  Fsuno  the  solar  flux  through  the  band  pass  of  the  sensor,  and  a  is  the  phase  angle. 
As  mentioned  previously,  light  curve  data  (Eq.  (lb))  is  sensitive  to  CdijjA  and  angles  data  (Eq.  (la))  is  sensitive  to 
CrA/m.  The  C,  coefficient  is  a  function  of  the  Cdijf  of  the  SO  surface.  In  general,  Cr  =  1  +  r  for  the  cannonball  model 
[14]  where  r  is  a  modeling  parameter  related  to  the  albedo.  Therefore,  to  allow  for  mass  to  be  separated  from  area, 
an  estimate  of  albedo  must  be  determined.  Note  that  r  and  Cdijf  are  not  the  same  and  rather  r  is  used  to  model  the 
effect  of  Cdijf  on  the  SRP.  It  is  shown  in  Ref.  [18]  that  C,  =  1  +  2Cdijj/3  for  a  diffuse  sphere  and  Cr=  1  +  4C dig/9  for  a 
diffuse  flat  plate  assuming  Lambert’s  cosine  law  with  the  normal  direction  aligned  with  the  Sun  direction. 

From  Eq.  (1)  and  the  relationship  between  C,  and  Cdijf,  it  is  clear  that  mass  separated  from  area  highly  depends  on 
the  assumptions  made  about  shape  and  surface  reflection.  Also  it  is  clear  that  brightness  and  the  SRP  force  are 
related  by  the  physics  that  they  model.  Therefore,  to  obtain  realistic  and  reliable  estimates  for  mass  both  the  SRP 
and  light  curve  models  must  be  physically  consistent.  Toward  this  end  a  general  SRP  model  is  defined  as 
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where  this  model  calculates  the  acceleration  due  to  N facet  facets.  The  terms  A*,  /a,  and  Fi(X)  are  the  area  of  the  kth 
facet,  fraction  of  facet  that  is  illuminated,  and  the  wavelength  dependent  solar  flux  respectively.  Additionally, 
(x)+  =  xH(x)  where  H(x)  is  the  Heaviside  step  function  which  is  unity  for  positive  values  and  zero  elsewhere.  The 
BRDF  (/,)  defines  how  light  is  reflected  from  an  opaque  surface  with  a  given  surface  normal  direction  ( N , )  based 

on  a  light  source  in  the  L  direction.  Finally  each  differential  dOjhp  spherical  element  of  outing  light  rays  has  a 
V  vector  defining  the  direction  of  the  acceleration  due  to  ray  exiting  through  this  spherical  element.  Then  the 
acceleration  can  be  calculated  by  defining  fr  depending  on  the  surface  properties  and  then  evaluating  the  integral  in 
Eq.(2). 

2.1  Shape  Model 

The  shape  model  considered  in  this  work  consists  of  a  finite  number  a  flat  facets.  For  curved  surfaces  this  model 
becomes  more  accurate  as  the  number  of  facets  is  increased.  Each  facet  has  a  set  of  three  basis  vectors 
(uf ,uf ,uB)  associated  with  it  as  defined  in  Figure  1(a).  The  unit  vector  uf  =N  points  in  the  direction  of  the 
outward  normal  to  the  facet,  and  the  vectors  uf  and  uf  are  in  the  plane  of  the  facet.  The  notation  superscript  B 
denotes  that  the  vector  is  expressed  in  body  coordinates.  The  SOs  are  assumed  to  be  rigid  bodies,  so  the  unit  vectors 
u  f ,  uf  and  uf  do  not  change  in  time  since  they  are  expressed  in  the  body  frame. 


The  light  curve  and  the  SRP  models  discussed  in  the  next  sections  require  that  these  vectors  be  expressed  in  inertial 
coordinates.  Since  the  SO  body  is  rotating,  these  vectors  will  change  inertially.  The  body  vectors  can  be  rotated  to 
the  inertial  frame  by  the  standard  attitude  mapping  given  by 

uf  =  ^(qf  K  k  =  u,v,n  (3) 

where  /\(q  ( )  is  the  attitude  matrix  mapping  the  inertial  frame  to  the  body  frame  using  the  quaternion 
parameterization.  The  superscript  I  denotes  that  a  vector  is  expressed  in  inertial  coordinates.  The  unit  vector 
ufiti  =  L  points  from  the  SO  to  the  Sun  direction,  and  the  unit  vector  u(bs  =  V  points  from  the  SO  to  the  observer 

as  shown  in  Figure  1(a).  The  vector  uf  =  His  the  normalized  half  vector  between  u(ln  and  u(bs  also  shown  in 
Figure  1(a).  This  vector  is  also  known  as  the  Sun-SO-Observer  bisector  and  has  an  angle  a  from  the  normal 
( cos  a  =  N .  ■  H  ).  Each  facet  has  an  area  Aa  associated  with  it.  Once  the  number  of  facets  has  been  defined  and 
their  basis  vectors  are  known,  the  areas  A*  define  the  size  and  shape  of  the  SO. 
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Figure  1.  Geometry  and  Shape  Description 

For  the  development  of  the  measured  light  curve  data,  faceted  SO  rectangular  cuboid  shape  models  as  shown  in 
Figure  1(b)  are  used.  The  rectangular  cuboid  model  is  described  by  three  parameters,  Z,  a ,  and  d ,  which  are  the 
length,  width  and  height,  respectively. 

2.2  Review  of  Bidirectional  Reflectance  Distribution  Functions 

A  number  of  bidirectional  reflectance  distribution  functions  (BRDFs)  exist  in  literature.  These  models  are  based  on 
the  BRDF  that  models  light  distribution  scattered  from  the  surface  due  to  the  incident  light.  The  BRDF  at  any  point 
on  the  surface  is  a  function  of  two  directions,  the  direction  from  which  the  light  source  originates  and  the  direction 
from  which  the  scattered  light  leaves  the  observed  the  surface.  These  models  include  Ashikhmin-Shirley  [19],  a 
simplified  Blinn-Phong  [20]  and  Cook-Torrance  [21],  The  BRDF  models  how  light  energy  reflects  off  of  surfaces 
and  how  this  reflected  energy  is  distributed. 

The  model  in  Ref.  [19]  decomposes  the  BRDF  into  a  specular  component  and  a  diffuse  component.  The  two  terms 
sum  to  give  the  total  BRDF: 


The  diffuse  component  represents  light  that  is  scattered  independent  of  the  angle  of  incident  light,  and  the  specular 
component  represents  light  that  is  concentrated  about  some  direction  (mirror-like).  Reference  [19]  develops  a  model 
for  continuous  arbitrary  surfaces  but  simplifies  for  flat  surfaces.  This  simplified  model  is  employed  in  this  work  as 
shape  models  are  considered  to  consist  of  a  finite  number  of  flat  facets.  Therefore,  the  total  observed  brightness  of 
an  object  becomes  the  sum  of  the  contribution  from  each  facet: 

fr={dRd  +  sRs)  (5) 


which  depends  on  the  diffuse  bidirectional  reflectance  ( R,i)  and  the  specular  bidirectional  reflectance  (Rs)  and  the 
fraction  of  each  to  the  total  ( d  and  .v  respectively  where  d  +  s  =  1).  These  bidirectional  reflectances  are  calculated 
differently  for  the  various  models.  In  each  model,  however,  p  is  the  diffuse  reflectance  (0  <P<  1),  and  Fo  is  the 
specular  reflectance  of  the  surface  at  normal  incidence  (0  <  Fo  <  1).  To  be  used  as  a  prediction  tool  for  brightness 
and  radiation  pressure  calculations,  an  important  aspect  of  the  BRDF  is  energy  conservation.  For  energy  to  be 
conserved,  the  integral  of  the  BRDF  times  cos  (0r)  over  all  solid  angles  in  the  hemisphere  with  0r  <  90  needs  to  be 
less  than  unity,  with 

2xX/2 

J  j  fr  cos(0,.)sin {er)derd(f>  =  Rd  +  Rs  (6) 

0  0 

For  the  BRDF  given  in  Eq.  (5),  this  corresponds  to  constant  values  of  Rc/  =  dp  and  R,  =  sFo.  The  remaining  energy 
not  reflected  by  the  surface  is  either  transmitted  or  absorbed.  In  this  paper  it  is  assumed  the  transmitted  energy  is 
zero.  A  review  of  the  various  BRDF  models  is  provided  in  Table  1. 


Ashikhmin-Shirley:  In  addition  to  cl  p,  and  Fo,  the  Ashikhmin-Shirley  BRDF  has  two  exponential  factors  (nu,  nv ) 
that  define  the  reflectance  properties  of  each  surface.  The  Ashikhmin-Shirley  diffuse  and  specular  reflectivities  are 
not  constant,  however,  but  rather  complicated  functions  of  illumination  angle,  exponential  factor,  and  the  diffuse  and 
specular  reflectances.  In  all  cases,  however,  Rd  +  Rs<  1,  so  energy  is  conserved. 

Blinn-Phong:  The  specular  bidirectional  reflectance  of  the  original  Phong  model  is  proportional  to  (Si  -r)  ,  where 

R  is  the  perfect  mirror-like  reflection  of  L  .  Blinn  [20]  proposes  that  H  be  used  instead  of  R  to  make  it  easier  and 
faster  to  calculate.  Unfortunately,  both  versions  of  the  model  do  not  conserve  energy  and  thus  are  unsuited  for  the 
purposes  of  brightness  estimation.  The  model  can  be  made  to  conserve  energy,  however,  by  modifying  the  leading 
term.  In  keeping  with  the  desire  for  simplicity  in  this  model,  the  leading  term  is  chosen  to  only  be  a  function  of  the 
exponential  factor  and  set  to  yield  a  reflectivity  equal  to  the  mirror-like  reflection  of  Eq.  (3)  at  normal  illumination. 
In  addition  to  d,  p,  and  Fo,  the  simplified  Blinn-Phong  BRDF  has  a  single  exponential  factor  (n)  that  defines  the 
reflectance  properties  of  each  surface. 

Cook-Torrance:  This  model  has  the  facet  slope  distribution  function  ( D ),  the  geometrical  attenuation  factor  (G)  and 
the  reflectance  of  a  perfectly  smooth  surface  ( F)  with  g  =  n2  +  ( V  •  H  )2  -  1  and  the  index  of  refraction 
n  =  (l  +  1  —  -J  l]  i  )•  I*1  addition  to  cl  p,  and  Fo,  the  Cook-Torrance  BRDF  has  a  facet  slope  (m)  parameter  that 

defines  the  reflectance  properties  of  each  surface.  The  facet  slope  parameter  of  the  Cook-Torrance  BRDF  and  the 
exponential  factor  of  the  Ashikhmin-Shirley  and  Blinn-Phong  BRDFs  are  roughly  related  by  n  =  2/m2. 


Table  1.  Review  of  Various  BRDF  Models 


3.  SOLAR  RADIATION  PRESSURE  DEPENDENCE  ON  BRDF 

As  described  in  Ref.  [22],  in  general,  R,i  and  Rs  are  complicated  functions  of  illumination  angle  and  material 
properties  represented  by  parameters  within  the  particular  BRDF  model.  For  certain  BRDFs,  however,  the  integral 


can  be  solved  analytically.  For  example,  the  BRDF  with  a  Lambertian  diffuse  component  and  purely  “mirror-like” 
specular  component  (where  R  =  2(l  ■  n)n  -  L  is  the  direction  of  mirror-like  reflection  of  L  )  is  given  by 
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where  8  is  the  Delta  function.  If  the  BRDF  of  Eq.  (7)  is  substituted  into  Eq.  (2),  the  following  acceleration  due  to 
the  SRP  is  found  by 
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For  a  more  complicated  BRDF,  the  exact  solution  obtained  by  numerically  integrating  Eq.  (2)  is  different  than  the 
idealized  solution  obtained  by  Eq.  (8).  The  numerical  integration  of  Eq.  (2)  over  the  full  hemisphere,  however,  is 
time  consuming  and  might  be  prohibitive  to  calculate  in  certain  applications.  In  this  paper,  correction  factors  for  Eq. 
(8)  from  Ref.  [22]  are  used  and  the  acceleration  is  calculated  in  a  single  step  by  using 
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where  the  A  factors  depend  on  the  illumination  angle  and  parameters  within  the  BRDF  model.  Thus,  for  A,;  =  1,  the 
diffuse  reflectance  is  Lambertian  and  for  Avi  =  Av2  =  1 ,  the  specular  reflectance  is  mirror-like.  The  further  these  A 
factors  are  from  unity,  the  greater  the  difference  between  the  BRDF  model  and  the  idealized  BRDF  of  Eq.  (8). 

4.  UNSCENTED  KALMAN  FILTER  FORMULATION 

Applying  the  UKF  structure  for  attitude  estimation  has  some  challenges.  For  instance,  although  three  parameter  sets 
are  attitude  minimal  representations,  they  inherently  have  singularities.  On  the  other  hand  the  quaternion 
representation,  which  is  a  four  parameter  set  with  no  singularity,  has  a  nonlinear  constraint  which  results  in  a  nearly 
singular  covariance  matrix.  This  does  not  allow  use  of  quaternions  in  a  straightforward  UKF  implementation. 
Reference  [23]  overcomes  these  challenges  by  utilizing  generalized  Rodrigues  parameters  (GRPs),  a  three  parameter 
set,  to  define  the  local  error  and  quaternions  to  define  the  global  attitude. 

Complete  explanations  of  the  quaternion  and  its  mapping  to  GRPs  are  provided  in  Refs.  [24]  and  [25].  In  the  UKF 
implementation  described  in  Ref.  [23],  the  covariance  matrix  is  interpreted  as  the  covariance  of  the  error  GRP 
because  for  small  angle  errors,  the  error  GRP  is  additive  and  the  UKF  structure  can  be  used  directly  to  compute 
sigma-points.  The  error  GRP  sigma  points  are  converted  to  error  quaternions  and  then  to  global  quaternions  for  the 
propagation  stage.  To  compute  the  propagated  covariance  the  global  quaternions  are  converted  to  error  quaternions 
and  then  back  to  error  GRPs.  The  process  is  then  as  follows:  error  GRP  — »  error  quaternion,  error  quaternion  — > 
global  quaternion,  global  quaternion  — *  error  quaternion,  and  finally  error  quaternion  — >  error  GRP. 

The  estimator  discussed  in  this  work  includes  constrained  states  since  the  area  states  are  constrained  to  be  positive. 
To  account  for  these  constraints  two  methods  are  used.  The  constrained  UKF  method  discussed  in  Ref.  [26]  is  used 
and  proxy  states  method  discussed  in  Ref.  [27]  is  used.  Reference  [26]  incorporates  the  information  of  the 
constraints  in  the  UKF  algorithm  by  projecting  sigma  points  and  estimates  onto  the  constraint  boundary  when  sigma 
points  or  estimates  are  found  outside  the  feasible  region.  Since  r  >  0  we  have  the  following  constraints: 

rA>  0  (10a) 


A  >  0 


(10b) 


Then  if  the  constraints  are  violated  during,  sampling  the  sigma  pointing  from  covariance  matrix,  after  propagation, 
or  after  the  update  step,  the  estimates  and  or  sigma  points  that  violate  the  constraints  are  projected  onto  the 
constraint  boundary. 

The  proxy  state  method  for  dealing  with  constraints  first  used  in  Ref.  [27]  is  based  on  defined  unconstrained  proxy 
values  and  estimating  those  values  instead  of  the  constrained  states.  To  account  for  the  constraints  within  the  filter, 
unconstrained  proxy  values  are  used  in  the  state  vector  and  estimation  filter  and  these  proxy  values  are  converted 
back  to  the  surface  parameter  value  when  needed.  Then  given  a  constrained  parameter  (ai  :  a  j  >  0,  a 2  :  0  <  a 2  <  1 ) 
proxies  bj  and  b:  can  be  defined.  The  conversion  equations  to  the  proxy  value  and  from  the  proxy  value  for  each  of 
the  surface  parameters  are 

bx  =  ln^q)  ax  =  exp(£>j )  (11a) 
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Depending  on  the  kind  of  state  constraint  either  Eq.  (1  la)  or  Eq.  (lib)  is  used. 

5.  SIMULATION  RESULTS 

Two  simulations  cases  are  presented:  an  Ashikhmin-Shirley  model  and  a  Cook-Torrance  model.  Both  cases  have 
the  same  dynamic  configuration,  with  same  initial  attitude,  true  orbit  and  true  angular  velocity.  The  cases  only 
differ  in  the  BRDF  model  used  to  simulate  and  process  the  measurements.  For  the  generation  of  data  for  the  true 
model,  an  equatorial  ground  station  is  chosen  as  the  site  of  the  observer.  Also,  all  scenarios  use  the  same  shape 
model  which  contains  six  sides  where  the  areas  and  BRDF  parameters  are  assumed  to  be  the  same  for  each  side. 
The  SO  is  simulated  to  fly  in  a  near-geosynchronous  orbit  in  a  trajectory  that  is  continuously  lighted.  This  is 
accomplished  by  inclining  the  orbit  by  30  degrees  and  choosing  an  appropriate  time  of  the  year,  thereby  avoiding  the 
shadow  cast  by  the  Earth.  The  simulation  parameters  are: 

•  Geographic  position  of  the  ground  site  is  0°  North,  172°  West  with  0  km  altitude 

•  Initial  inertial  position  and  velocity  are  chosen  as  r/=  [-7.893 lx  102  3.6679xl04  2.1184xl04]T  km 

•  Initial  inertial  velocity  is  chosen  as  v'  =  [-3.0669  -4.9425xl0“2  -2.8545xl0_2]T  km/s 

•  The  initial  time  of  the  simulation  is  May  8,  2007  at  5:27:55 

•  Initial  quaternion:  =  1 1/V2  0  0  1/V2]T 

•  A  constant  rotation  rate,  defined  as  the  body  rate  with  respect  to  the  inertial  frame,  represented  in  body 
coordinates,  is  used  given  by  q®  =  [0  0.00262  0]T  rad/s 

For  all  simulations  scenarios,  measurements  are  produced  using  zero-mean  white-noise  error  processes  with 
standard  deviation  of  0.1  for  magnitude  and  standard  deviation  of  1  arc-seconds  for  azimuth  and  elevation.  The 
initial  errors  for  the  states  are  10-30  deg  for  all  three  attitudes,  14  deg/hr  for  the  rotational  rates,  100  km  and  1  km/s 
for  the  position  and  the  velocity  errors,  respectively,  600  kg  for  mass,  and  0.1  for  both  the  diffuse  reflectance  and 
specular  reflectance.  The  initial  condition  error-covariance  values  are  set  to  (6.67)2  deg2  for  all  three  attitude 
components,  242  (deg/hr)2  for  the  rotational  rates,  and  1002  km2  and  l2  (km/s)2  for  the  position  and  the  velocity 
errors,  respectively.  The  initial  condition  error-covariance  values  for  the  mass,  diffuse  reflectance  and  specular 
reflectance  are  3002  kg2,  0.22,  and  0.22,  respectively.  The  time  interval  between  the  measurements  is  set  to  30 
seconds.  Data  is  simulated  for  20  nights  (about  20  orbits)  where  observations  of  the  SO  are  made  each  night  for  1 
hour.  The  simulation  results  are  plotted  versus  number  of  data  samples  since  there  are  large  time  gaps  between  each 
1  hour  data  arc. 

For  the  development  of  the  measured  light  curve  data  a  six  sided  facet  SO  shape  model  is  used.  The  three 
parameters  for  the  shape  model  are  given,  /,  a,  and  d  which  are  the  length,  width,  and  height  respectively,  as  shown 
in  Figure  1(b).  For  the  truth  model  a  six-facet  rectangular  shape  model  is  used  with  dimension  l  =  4  m,  a  =  2  m  and 
d  =  8  m.  The  mass  for  the  true  SO  is  selected  to  be  m  =  600  kg,  which  gives  a  maximum  area-to-mass  ratio  of  0.09 
for  this  model. 


(a)  Mass  Error  Estimates 


(b)  p  Estimate 


Figure  2.  Ashikhmin-Shirley  Results 

The  Ashikhmin-Shirley  and  the  Cook-Torrance  model  mass  estimation  results  are  shown  Figs.  2  and  3.  The  mass 
estimates  along  with  the  surface  parameter  estimates  are  also  shown  in  Figs.  2  and  3.  From  these  figures  it  can  be 
seen  that  the  surface  parameters  are  estimated  well  for  both  models  but  the  parameter  p  does  not  show  good 
observability.  This  is  due  to  the  fact  that  truth  surface  parameters  used  are  not  very  specular  and  therefore  the 
measurements  do  not  provide  good  observability  with  respect  to  p.  From  these  figures  it  can  also  be  seen  that  the 
Cook-Torrance  model  mass  estimates  coverages  faster  than  that  of  the  Ashikhmin-Shirley  model.  This  can  be 
attributed  to  the  fact  that  the  SRP  acceleration  calculated  from  the  Cook-Torrance  is  larger  then  that  of  the 
Ashikhmin-Shirley  model. 

6.  CONCLUSIONS 

In  this  paper,  an  unscented  Kalman  filter  estimation  scheme  using  light  curve  and  angles  data  was  presented  and  was 
used  to  estimate  mass  and  bidirectional  reflectance  distribution  function  (BRDF)  parameters  of  a  space  object  (SO) 
along  with  its  associated  rotational  and  translational  states.  This  work  extends  previous  work  by  using  physical 
consistent  solar  radiation  pressure  models  based  on  BRDFs.  These  models  allow  for  fusion  of  light  curve  and  angles 
data  where  the  parameters  used  in  both  models  are  consistent  and  therefore  provide  realistic  estimates.  This  work 
uses  an  assumed  shape  model  of  six  sides  and  estimated  the  area  and  BRDF  parameters.  Both  the  Ashikhmin- 
Shirley  and  the  Cook-Torrance  BRDF  models  were  used  to  simulate  and  process  the  light  curve  and  angle 
measurements.  Each  model  was  used  to  estimate  mass  along  with  the  other  aforementioned  states,  and  good 
performance  was  shown.  Using  an  unscented  filter  to  employ  brightness  magnitude  and  angles  data,  the  estimator 
was  able  to  determine  the  mass  of  a  SO  to  within  high  accuracy. 


(a)  Mass  Error  Estimates 


(b)  p  Estimate 
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Figure  3.  Cook-Torrance  Results 
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